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Abstract. The global dynamics of the ionized and neutral components in the 
environment of Io plays an important role in the interaction of Jupiter’s corotating 
magnetospheric plasma with Io. The stationary simulation of this problem was done in 
the MHD and the electrodynamics approaches. One of the main significant results 
from the simplified two-fluid model simulations (Saur et al., 2002) was a 
production of the structure of the double-peak in the magnetic field signature 
of the 10 flyby that could not be explained by standard MHD models. In this 
paper, we develop a method of kinetic ion simulation. This method employs 
the fluid description for electrons and neutrals whereas for ions multilevel, drift-kinetic 
and particle, approaches are used. We also take into account charge-exchange 
and photoionization processes. Our model provides much more accurate 
description for ion dynamics and allows us to take into account the realistic 
anisotropic ion distribution that cannot be done in fluid simulations. The first 
results of such simulation of the dynamics of ions in the Io’s environment are discussed 
in this paper. 
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1. Introduction 

The global interaction of the Jovian plasma torus with the Io is a fundamental 
problem in magnetospheric physics. It requires the solution of a highly nonlinear 
coupled set of integro-MHD/kinetic-Boltzmann equations which describe the dynamics 
of Jupiter’s corotating magnetospheric plasma, pickup ions and the ionosphere together 
with the atoms from Io’s environment. To leading order, the plasma and neutral atoms 
are coupled by resonant charge exchange, although other coupling processes are present. 
The characteristic scale of the ionized components is determined usually by the typical 
ion gyroradius, which for Io is much less than characteristic global magnetospheric scales 
of interest. By contrast, the mean free path of neutral particles and the typical ion 
gyroradius are comparable to characteristic exospheric scales such as the thickness of 
the exosphere and atmosphere, the distances separating the possible magnetic barriers 
and the surface of Io, etc. Consequently, the Knudsen number Kn = \/L exo (A, the 
mean free path of neutral particles and L exo , a characteristic exospheric scale), which is 
a measure of the distribution relaxation distance, satisfies Kn as 1. Thus it is difficult 
to assume that the distribution of the neutral atoms and ions can relax to a Maxwellian 
distribution, and one needs ideally to solve a Boltzmann equation for the neutral 
and ionized components in which charge exchange and photoionization processes are 
included. 

Several approaches for including the neutral component and pickup ions self- 
consistently in models describing the interaction of plasma torus with Io have been 
formulated. Early theoretical work was often done either in the context of a “thin” 
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atmosphere (e.g., see Cloutier et al., 1978) indicative of the surface temperature ( 130 K), 
or “thick” extended neutral atmosphere (e.g., see Goertz, 1980) more indicative of 
volcanic temperatures ( 1000K). Subsequent evidence (see the review by Lellouch, 1996) 
seems to indicate a mixed picture of the global atmosphere, which has a large extended 
corona like a thick atmosphere, but appears to be dominated by local major injection of 
hot (high speed) gas/dust plumes to high altitudes but only near active volcanic vents. 
Therefore, although the atmosphere is probably only locally thick, it still has a large 
extended neutral corona which might provide a sufficient source of impact ionization 
and photoionization to explain the plasma torus. 

Southwood et al. (1980) examined data from several Voyager instruments and 
examined the possible role of an intrinsic magnetic field for Io as a way to retain a 
robust enough ionosphere, which could provide enough conductivity for completing the 
Io-Jupiter current circuit. Neubauer (1980) presented an analytical model of 
the Alfven standing wave current system which connects current through 
the ionosphere of Io. Southwood and Dunlop (1984) and Ip (1990) suggested models 
stating that mass loading effects should result in the formation of a tail-like structure in 
the wake behind Io; as a consequence of the enhanced plasma density, the magnetic field 
perturbations are continued into the wake. Thus the field-aligned current is not only 
generated by Io itself, but also in the wake. Several years after Voyager, 3D numerical 
studies of the plasma flow past Io were performed using electrodynamic (Wolf-Gladrow 
et al., 1987), magnetohydrodynamic (Linker et al. 1989; Linker et al. 1991) and resistive 
magnetohydrodynamic (Kopp, 1996; Kopp, Birk and Otto, 1998) approaches. 
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There have been recent efforts to improve and extend the pre-Galileo simulations 
both in terms of the MHD (Combi et al., 1998; Linker et al., 1998; Kabin et ah, 2001) 
and the electrodynamic (Saur et al., 1999) approaches. These two approaches are 
distinguished by the physical assumptions which they each do and do not (or in some 
cases, can and cannot) include. However, MHD cannot, at least yet, include the effects 
of realistic conductivities (Hall and Pederson) or charge separation effects which are 
likely to be important very close to Io where the neutral densities are large and electric 
potential can introduce non-symmetric flow around body. They either include constant 
artificial conductivity (Linker et., 1998) or assume perfect conductivity (Combi et al., 
1998), however comparisons of the sets of published results do not indicate that this 
choice has any important consequence. 

The non-magnetic models produce magnetic field perturbations that are similar to 
the Galileo measurements, but none are quite as deep or as broad, and none have the 
reversal of the perturbation (the double-peaked structure) in the center of the wake. The 
magnetized models of Linker et al. (1998) produced a broad and deep perturbation, but 
not the self-reversal at the center of the wake (the double-peak or bite-out). Thus, the 
observational data, e.g., density and magnetic field profile along the Galileo trajectory 
cannot be fully explained by MHD models. But the most recent evidence from all the 
Galileo flyby data is that Io does not possess any substantial internal magnetic field of 
its own (Kivelson et al., 2001). 

In the electrodynamical model of Saur et al. (2002), they added an extra ionization 
source, based on the high-energy bidirectional electrons observed by Williams et al. 
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(1996,1999) and Frank and Paterson (1999). These high-energy electrons were included 
as an energy source in addition to electron impact ionization by the thermal electrons 
and photoionization to create a dense plasma wake as is observed. Simulations without 
these extra bidirectional electrons show a nearly empty wake Saur et al. (1999). They 
also concluded that the electrodynamic part of Io’s interaction is best described as an 
ionosphere-like interaction rather than a comet-like interaction (Saur et al., 2003). Saur 
et al. (2002) demonstrated first that the diamagnetic and inertia currents 
are responsible for formation of the double-peak magnetic profile along the 
10 trajectory. They received the magnetic field profile that has an oscillating 
structure and the maximum in the magnetic field profile is to some extend 
was a narrower than in observation for standard model of the atmosphere 
(Fig. 11, Saur et al., 2002). In a case of longitudinally symmetric atmosphere 
they received a strong double-pick structure but the spatial scale of this 
structure was a twice smaller than in observation (Fig. 14, Saur et al., 2002) 
and the space scale of these peaks are much smaller than in observation. 

In this paper, we describe a new approach to solving the time-dependent Boltzmann 
equation together with a hybrid plasma (ion kinetic) model in three spatial dimensions 
using a prescribed atmosphere model for Io. A Boltzmann simulation is applied to model 
charge exchange between (incoming and pickup) ions and immobile atmosphere. Several 
simulations are run, and the results described. We show for the first time the predicted 
distribution for ions and the electromagnetic field throughout Io’s environment. The 
results of such kinetic simulations are compared with those obtained from a related 
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MHD model and observational data. We also found the best parameters for our model 
that describe satisfactory the main behavior of the observational data at least on the 
qualitative level. 

2. Formulation of the Problem and Mathematical Model 

2.1. Simulation Model 

To study the interaction of the plasma torus with the ionized and neutral 
components of Io’s environment we use quasineutral hybrid models for ions and 
electrons. The interaction of the ions with atmosphere is dominated by charge exchange. 
The atmosphere is considered to be an immobile component in this paper. The general 
scheme of the global interaction of the plasma torus with Io and the Galileo 10 trajectory 
is given in Fig. 1. The Galileo 10 flyby occurred nearly in the equator plane of Io and 
perpendicular to the direction of the plasma wake defined by the corotating plasma flow 
past Io. The spacecraft trajectory passed approximately 900 km down-stream of Io (in 
the sense of the plasma torus flow). In our coordinate system the Z axis is parallel to 
Uq (corotational plasma velocity), Y is aligned with the spin axis, and X = Y x Z. The 
relation between our coordinate system and IphiO (X*,Y*,Z*) system (see Fig. 3 from 
Kivelson et al., 2001) is the following: X is parallel to Y*, Y is aligned with Z*, and Z 
is parallel to X* . 

In the hybrid simulations described here, the dynamics of upstream ions and 
implanted ions is described in a kinetic approach, while the dynamics of the electrons is 
described in a hydrodynamical approximation. 
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The single particle ion distribution function / s (t,x,v) has to fulfill the Vlasov 
equation 

d , d , F5, 

where F symbolizes all forces acting on the ions. 

The single ion particle motion is described by the equations 


d r 


S,l 


= v S) b 


dv 3i , _ Z a e / v 3 j x B\ 
dt M s V + c ) ’ 


E* = E - 


( 2 ) 


dt dt 

Here and M„ denote the charge state and the mass of the ions, the subscript 
s denotes the ion species (s = 1 for incoming ions and s = 2 for pickup ions) and the 
index l is the particle index. a e ff is an effective conductivity that may include Coulomb 
collisions and collisions due to particle-wave interaction, in the sense of a standard 
resistive Ohm’s law. Note that the resistivity used in Eq. (2) must depend on individual 
velocities of ions and electrons. However, we use the effective resistivity and an effective 
electric field E* in Eq. (2) to satisfy momentum conservation for the plasma system 
(Zueva et al., 1975). 

In the nonradiative limit Ampere’s law is given by 

47T 


V x B; 


( 3 ) 


and the induction equation (Faraday’s law) by 

15B 


c dt 


+ V x E = 0. 


( 4 ) 


The total current is given by 


J = J e + Ji; Ji = 53 Z s en s \J s , 

S = 1 


( 5 ) 
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here U s is the bulk velocity of ions of the type s. 

We further assume quasi-neutrality 

2 

n e = Yl Z s n s- ( 6 ) 

5=1 

For massless electrons the equation of motion of the electron fluid takes the form of 
standard generalized Ohm’s law (e.g. Braginskii, 1965): 

E = — (J e X B) - — Vp e + —Re, (7) 

en e c en e en e 

where p e = nm e (v ^)/ 3 = n e k B T e , and v' e are the scalar electron pressure and the thermal 
velocity of electrons, and the electron current is estimated from Eq. ( 5). 

The term including R e , symbolizes the mean momentum change of the electrons 
due to their collisions with the ions, which can be parametrized in terms of the collision 
frequency u es between electrons and ions of the species s as 

Re — -n e m e u es (U e - U s ) (8) 

Since we suppose that electron heating due to collisions with ions is very small, the 
electron fluid is considered adiabatic: 

Pe oc n 5 J 3 . (9) 

The ion kinetic approach allows us to take into account the effects of anisotropy of ion 
pressure, the correct mass loading processes, a penetration of ions across the ionosphere, 
and the asymmetry of plasma flow around the Io. Remember that the fluid models 
which account only for the scalar ion pressure may result in an extra- effusion of the 


pickup ions along the Alfven wing. 
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Charge exchange may be included in a hybrid model by a simple procedure (Lipatov 
et al., 1998; Lipatov, 2002). The total loss rate in s _1 for the ionized component has the 
form 

&x(x, V, t) = J / 0 (x, V tt , t)V re l, a a e x{Vrel,a)d 3 V a , (10) 

where the velocity of ions relative to an atom with velocity v a is V re ^ a — |v — v a |. 

If we suppose that the neutral component has a Maxwellian distribution, then {3 ex 
may be approximated as (Ripkin and Fahr, 1983) 

0ex(x, v, t) ~ n 0 (x, t)Vr e l, a <J ex (Vrel,a)- (11) 

In this paper we assume that the neutral component of Io’s atmosphere has zero bulk 
and thermal velocities. If one needs to include the nonzero temperature and bulk 
velocity, the effective average velocity of ions relative to an atom with velocity may be 
found from the general equation of Ripkin and Fahr (1983). 

Let the time interval t* with respect to the charge exchange be a random variable 
with a distribution function 

wi(t*) = exp[— f Pexjdt], (12) 

Jto 

where / is the particle index. 

A survival probability against the charge exchange event w ex may be written as 
follows: 

Wex,l(t) = ex P[~ / Pex-Jtdt]. (13) 

Jto 

The integration is over the trajectory of the particle with index l. At the time of 
creation (either at the boundary of the calculation box, or at the moment of charge 
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exchange), an ionized particle has initial coordinates x;(t 0 ) = x ij0 , V((t 0 ) = v; ]0 , a weight 
a t (to ) = aifl, and a survival probability w ex j(t Q ) = 1. For each new ion with index l, we 
have to determine the critical probability w* ex l when charge exchange will occur, and 
this is done using the relation 

<,/ = 6 ( 14 ) 


where £ is random number on the interval from 0 to 1. During the calculation we have 
to identify those particles for which the probability of survival satisfies the condition 


^ex,i < W* ex l . (15) 

If the particle satisfies the condition (15), then we have to exchange the velocity of this 
ion with the velocity of an atom from the atmosphere of Io. In the present simulations, 
we do not take the cross section a into account in (10) because of the weak dependence 

on \v — v p \ (as was done in the work of Malama (1991)). If charge exchange occurs, then 

a new neutral particle begins its motion with v* = 0 and w eXf t = 1. 

Io’s environment model includes 2 sources of pickup ions: the extended coma halo 
distribution with pickup ion production law 

G OC ^i^atmosbke^r (lb) 


for 1 x ri 0 < r < 7 x r Io and the exosphere distribution close to Io with ion production 
law 

G oc t'in atmos W int exp[— (r - r Xo ) / H^ mos }. (17) 


Here W ext and W int denote the fraction of the total ion production rate in the halo 
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and exosphere (W tI , t — (100 — 95)% of the total), respectively. n atmos denotes the 
maximum atmosphere density, is effective ionization rate. 

The inner region of the Io’s ionosphere is described by immobile ions with the 
following density distributions: 

0.5n iono exp(— (r - r lo )/H efS ), for r > r Io , 

rciono(l - 0.5exp(-|r - r Io |/// e //)), for r < r io , 

We assume that the incoming flow has a small finite resistivity to suppress the “shot” 
noise fluctuations. We also take into account the effect of finite conductivity of Io’s 
body so that 

<7e// = CTup, for r > rio, 

°eff = CTio, for r < rio, 

Our code solves equations (1) - (8), (5) - (9), (10) - (15) and (16) - (17). 

Initially the computational domain contains only supersonic plasma torus flow 
with a homogeneous spatial distribution and a Maxwellian velocity distribution; the 
pickup ions have a weak density and spherical spatial distribution. The magnetic 
field and electric fields are B = Bo and E = Uo x Bo. Inside the Io electromagnetic 
fields are E =0 and B = Bo- In the cases examined in Figs. 9-16 we choose 
B 0 /Bo = (0,1,0); Hjo/Eq = (1,0,0) and in the cases examined in Figs. 2-8, 17-19 we 
choose B 0 /B 0 = (0.0394,-0.9854,0.1657) and E 0 /E 0 = (-0.98,-0.0394,0.0). 

At t > 0 we begin to inject the pickup ions with a distribution according to Eqs. 
(16-17). Far upstream (z = —DZ/ 2), the ion flux is assumed to have a Maxwellian 
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: distribution, 

i / = noo(7ru t 2 h )" 3/2 exp -- V ^ , (18) 

lV th . 

where v t h and U are the thermal and the bulk velocities of the plasma torus flow, 

respectively. 

i 

Far downstream, we adopt a free escape condition for particles and Sommerfeld’s 
radiation condition for the magnetic field. On the side boundaries (y = ±DY/2 and 

| 

x = ±DX/ 2), periodic boundary conditions are imposed for incoming flow particles and 
the electromagnetic field. In some cases we also tested the use of the upstream boundary 
condition for electromagnetic field on the side boundaries. In these situations we also 
employ a buffer zone with thickness about of 10 x Ax where a smoothing procedure 
provides a transition for electromagnetic fields from the perturbed value to the upstream 
value on the side boundaries (see e.g. Umeda et al., 2001), effectively allowing Io- 
generated Alfven disturbances to propagate away. The pickup ions come out from the 
computational domain when they intersect the surfaces x = 5 x Ax, x — DX — 5 x Ax, 
y = 5 x A y, y = DY — 5 x Ay, z = 5 x A z, z — DZ — 5 x A z. Thus the flux of 
the pickup ions is absent on the side boundaries. At Io’s surface the particles may be 
reflected or absorbed. Note that the position of the Io is x = 0, y = 0, z = 0. 

The three-dimensional computational domain has the dimensions DX = 20 L, 

DY = 20L, and DZ — 10L, or DX — 14 L, DY = 14 L, and DZ = 12 L, or DX = 10 L, 
DY = 10 L, and DZ = 8 L, where L equals the radius of Io, r\ Q = 1800 km. We use the 
meshes of 201 x 201 x 101, or 161 x 161 x 121, or 141 x 141 x 121 grid points, and 2 x 10 8 
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and 5 x 10 7 particles for ions and pickup ions, respectively, for a homogeneous mesh 
computation. The time step At satisfies the condition v max At < min(Ax, Ay, Az)/8. 

The relationship between dimensional (U, E , B , p e , n, T) and dimensionless ({/', 
E ! , B f , Pg, n 7 , T f ) parameters may be expressed via dimensional upstream values as 
follows: 

U = UVo, E = E'B 0 Uo/c, B = B'Bo, Pe = Pe'Pe o, 

n = n'n 0 , T = T'MJJl , (19) 

whereas the dimensional time and distance may be expressed via the bulk velocity f7o 
and characteristic scale L: 

t = t'L/Uo , ^ (20) 


2.2. Numerical Method 

We employ a standard particle-in-cell (PIC) method in the case of a homogeneous 
grid. The time integration of the particle motion equations uses a leapfrog scheme. The 
time integration of the electromagnetic field equations uses an implicit finite difference 
scheme (see, e.g., Lipatov (2002)). We use different time steps for particle and field 
pushing (subcycling). This code was optimized for parallel computation using MPI and 
OMP. 

Since the gyroradii must be resolved, a grid point spacing of less than 1 gyroradius 
is required in order to avoid numerical dispersion and dissipation. On the other hand, 
good statistics are required, therefore a sufficiently large number of particles per cell 
have to be used (i.e. low “shot” noise). We use a homogeneous mesh for the simulations 
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presented here. 

3. Results of the Simulation 

To study the interaction of the plasma torus with the ionosphere of Io the following 
set of Jovian plasma torus and ionosphere parameters were adopted in accordance 
with observational data (Frank et al., 1996; Kivelson et al., 1996): corotation velocity, 
Uq = 56.8 km/s; plasma density, no = 3500 cm -3 ; plasma temperature, 92 eV; mean ion 
mass, Mi = 22 M p ; ratio of specific heats, 7 = 5/3; ion and electron betas, $ = 0.039; 

= 0.0022; magnetic field, B = 1800 nT; Alfven and sonic Mach numbers Ma = 0.4 
and Ms = 2.2. 

For pickup ions we use the following parameters: total ion production rate, 

Q ion = (0.7 — 2.) x 10 28 s -1 ; mean ion mass, Mpi = 22 M p ; electron exosphere and 
ionosphere betas, /3 e ,ex 0 = 0 — 0.0022, /3 e ,iono = 0 - 0.0005; effective cross section for 
charge exchange, a exc h = 1-5 x 10 -15 cm 2 ; effective ionization rate, 17 = 10 -5 s -1 ; 
atmosphere scale height, H^ mos = 0.06 — 0.09 x ri 0 = 108 km— 162 km; maximum value 
of the density of the atmosphere, n atmos = (0.5 — 10) x 10 8 cm -3 ; effective ionosphere 
scale height, H e ff = (0.0001 — 0.00025) x r\ Q =(0.18-0.45) km. Note that in calculation 
the effective ionosphere scale height is smoothed over the nearest grid cell and its 
specific values is not crucial. However, for charge exchange between the ions and atoms 
we use the analytical formula for the density of the atmosphere without any smoothing. 
The effective dimensionless diffusion length of the upstream and ionosphere plasmas 
is ld,up = 0.00125, whereas for Io’s body the diffusion length is = 0.0125, where 


Id = 1 /Re and the magnetic Reynolds number is Re = 4nU Q Lo e f f/c 2 . The average 
Lunquist number for the ionosphere may be estimated as 
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Re 


M A 


where the characteristic diffusion and Alfven times are = AnoL 2 /c 2 and ta = L/v a- 
Note that Linker et al. (1998) used the following Lundquist number, Lq = 500 — 1000, 
magnetic Reynolds number, Re = 200 — 400 and diffusion length, l ( j = 0.00125 for 
the background plasma. This value of the diffusion length corresponds to an effective 
conductivity of about 1.2 x 10 8 f2 -1 m -1 far in upstream. Since many plasma and 
atmosphere parameters are still uncertain we have to study a wide spectrum of simple 
models in order to choose the best one for interpretation of the observational data. We 
can present here only a sample of the wide spectrum simulation results which are in 
agreement or in disagreement with observations in order to illustrate the dependence of 
the plasma environment near Io on the input parameters of the problem. The global 
structure of Io’s environment is determined by a set of dimensionless independent 
parameters such as M\, /3 P , /3 e , Mp\/M p , ion production and charge exchange rates, 
diffusion lengths, and the ion gyroradius e = p C i/L. For real values of the magnetic field 
the value of the ion gyroradius is about 8 km which is calculated by the use of the 
local bulk velocity. The dimensionless ion gyroradius and grid spacing have the values 
e = 0.0045 and A X /L = 0.1. In order to study the ion kinetic effect we have to resolve 
the ion gyroradius on the grid. For this reason we use the artificially increased value of 
this parameter, e = 0.126, that is about of a value of A X /L now. Although, such way 
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allows us to study the some kinetic effects for realistic physics we have to extrapolate 
our results into a realistic scale. By scaling the gyroradius in such a way, we preserve 
the ratio of M A , Ms, /3 p , /3 e , and more importantly accurately preserve any anisotropy 
of the ion distribution function with respect to the magnetic field. 

3.1. Global Structure of the Io’s Environment 

Let us consider first the global picture of the interaction of the plasma torus with Io 
in a case with an ionization rate <3j on = 3.03 x 10 27 s -1 (run ceie2, Table 1) and charge 
exchange of form ^ + B exp . The maximum value of the atmosphere density 
is n atm os = 5 x 10 7 cm~ 3 and diffusion lengths are Zd.up = 0.00125 and Z^io = 0.00125. 
Figures 2 and 3 demonstrate the asymmetrical distribution of the torus plasma ion (top) 
and pickup ion density (bottom) in the x-z and y-z planes. One can see the increase 
of the incoming ion density upstream of Io. The pickup ion motion is determined 
mainly by the electromagnetic drift. The motion along the magnetic field is due to 
the thermal velocity and the gradient of the electron pressure. The asymmetrical 
distribution of the incoming ions in the y-z plane may be explained by an existence 
of the B z component of the upstream magnetic field. The inclination of the magnetic 
field results in the asymmetrical boundary condition for ion dynamics (penetration and 
reflection) in the Io’s ionosphere and the asymmetrical Alfven wing. The incoming 
ions flow around the Alfven wing so that only a small portion of the incoming flux 
penetrates the wing, resulting in a decrease in the incoming ion density. This effect is 
stronger in upper half-plane (y > 0) (Fig. 3, top) because the Alfven wing has a stronger 
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front transition in this region than in a lower half-plane. The pickup ion distribution, 
Fig. 3 (bottom), gives the correct value for the inclination of the Alfven wing, 21-21.5 
degrees. The density profile is a little bit disturbed near the side boundaries, however, 
this perturbation does not affect the region close to Io as discussed already. Figures 4 
and 5 show the distribution of the electric and magnetic fields. The asymmetry of the 
distributions in E and B appears to be caused by finite gyroradius effects of incoming 
and pickup ions. A weak perturbation of the magnetic field was observed near the 
ionosphere of Io: compression of the magnetic field upstream and decompression in the 
plasma wake. 

Figure 5 also shows the formation of a strong Alfven (and whistler) wing in the 
direction of the magnetic field. The perturbation of the electric field inside the wings is 
very strong and it may affect ion dynamics so that particles flow around the wings. The 
formation of the Alfven wing in a subalfvenic flow near Io was studied first analytically 
by Neubauer (1980). A excitation of a whistler wave near a plasma cloud was studied by 
using 3-D hybrid simulation in Lipatov (2002). The above wave propagation is closely 
connected with the generation of low frequency waves by the harmonic dipole (local 
source) in the magnetized plasma. The first analytical studies of these effects may be 
found, for example in (Van’yan and Lipatov, 1972, and references therein). 

Figures 6 and 7 show the velocity arrows of incoming and pickup ion velocities. The 
incoming ions flow around the effective obstacle that is produced by pickup ions and 
the ionosphere. The pickup ions flow from the “corona” across the magnetic field due to 
electromagnetic drift whereas the motion along the magnetic field is determined by the 



20 


thermal velocity of ions and the electron pressure. Figure 8 shows the two-dimensional 
cross section for total ion density in the x-y plane. One can see the asymmetry of 
distribution relative to the x-axis due to the angle between the bulk velocity and the 
magnetic field upstream and to the y-axis due to effects of the finite ion gyroradius. 

3.2. Effect of Electron Temperature on the Plasma Environment 

Our model takes into account the three characteristic electron temperatures that 
do not play any role in past MHD models: a) the temperature of electrons in the plasma 
torus, Te (Up ; b) the temperature of electrons that are created together with pickup ions, 
T e pi; c) the temperature of electrons that are created with ions in the ionosphere, 

T e< iono- Note that in all cases discussed here the pickup ions are generated only in the 
exosphere, i.e W ext = 0. 

Let us consider the case (a). Figure 9 shows the one-dimensional cuts of the total 
ion density, temperature and magnetic field along the x-axis for z = 1.5 x ri 0 and 
y — 0. One can see two peaks in the distribution of the density, each of which has a 
width of about ri 0 and maximum densities of about 4 and 6.5 relative to the upstream 
density. The depletion of the density at x = 2 L may be explained by the following. 
When the supersonic flow passes around the obstacle, a wake with a decreased density 
is formed. The simple example of a such void is presented by the lunar plasma wake. 
However, in our case an asymmetrical pickup ion high density obstacle provides an 
asymmetrical void in Io’s plasma wake. The temperature profile has two peaks with 
maximum temperatures of about 2.3 and 0.7. These peaks do not correlate with the 
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peak in density because the temperature of the ions is determined by the heated ions 
from the incoming flow and pickup ions. The magnetic field profile shows a decrease 
in the magnetic field with a minimum value of about 0.5 in the plasma wake and with 
some significant oscillations (Fig. 9). Figure 10 shows a two-dimensional cross section of 
the total ion density across the plasma wake ( z = 1.5 x ri 0 ). One can see three strong 
maxima. This distribution is determined by finite gyroradius effects of ions. 

If we take into account the temperature of electrons that are created with pickup 
ions (case (b)) the distribution of the ions in the plasma wake may be changed 
significantly. The one-dimensional total density profile has one peak with a value of 
about 8 (Fig. 11). The temperature profile has also one strong peak with a value of 
about 1.9 and the magnetic field has a depletion with a minimum value 0.7 (Fig. 11). 
Figure 12 shows a two-dimensional cross section of the total density across the plasma 
wake. One can see a strong peak narrow in x-direction and wide in y-direction. 

In case (c), when the temperature of electrons inside the ionosphere is also taken 
into account, we have the following distribution of plasma parameters in the wake. 
Figure 13 shows one-dimensional cuts of the total density, temperature and magnetic 
field. The density profile has a wide peak with a thickness of about 4ri 0 and a peak 
value of n rs 6. The temperature profile has a two-peak distribution. The value of 
the temperature at the two peaks is about 0.6 and 0.55. The peaks are separated by 
a distance of 4 x ri 0 . The magnetic field profile has a minimum at the plasma wake 
with B = 0.55 (Fig. 13). Figure 14 shows the two-dimensional cross section of the 
total density across the plasma wake. One can see a strong peak which is narrow in 
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the y-direction and wide in the ^-direction. By introducing cool electrons inside the 
ionosphere we are trying to account in an appropriate way for the cooling effect on 
electrons through electron-collisions close to Io. 

Finally, we show the two-dimensional cuts of pickup ion density for all three cases, 
Fig. 15 (top, middle, bottom). One can see that the temperature of electrons which 
are created in different regions of Io’s environment may affect strongly the pickup ion 
distribution. Figure 15 demonstrates the importance of the electron pressure in the 
region close to the Io. 

In the case of the absence of pressure of pickup and ionospheric electrons, the 
interaction of the plasma torus ions with Io is strongly asymmetrical in the x-z plane 
due to the finite ion gyroradius. In the plasma wake a ’’space-mixing” of the pickup ions 
is observed that results in the formation of an asymmetrical tongue-type distribution 
in the x-z plane. The effusion of pickup ions in the y direction is also weak because of 
small values of A p e and the velocity of the pickup ions in y direction, Fig. 15, (top). 

In the opposite case, when the gradient of the pressure of the pickup electrons and 
the pressure of ionospheric electrons cause a strong electrostatic field, the pickup ion 
density distribution in the x-z plane becomes more symmetrical and wider across the 
wake compared with the above case (cf. Fig. 15 (top) and (bottom)). The gradient of 
the electron pressure also causes a strong effusion of the pickup ions along the Alfven 
wing (see, Fig. 15 (bottom) and (top)). 

Finally, in the intermediate case when the pressure of ionospheric electrons is small 
the asymmetrical distribution of the pickup ion density with a strong effusion along the 
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Alfven wing is observed, Fig. 15 (middle). So, we can conclude that the pressure of the 
pickup ions is responsible for effusion processes while the pressure of the ionospheric 
electrons is responsible for the formation of the symmetrical distribution with the 
maximum in density located near the (y=0) plane. 

All three cases, (a), (b) and (c), demonstrate strong asymmetries in the distributions 
of pickup ion density, temperature and magnetic field in comparison with MHD and 
electrodynamic models. 

3.3. Effect of Pickup Ion Injection Distribution on the Plasma Environment 

To study the role of pickup ions from the halo we simulated Io’s environment 
for cases when the ion production rate in the halo corresponds to W ext — 10% and 
W ext = 30% of the total (halo plus exosphere) ion production rate. Note that all other 
parameters are the same in these cases. 

Figure 16 (solid line) presents one-dimensional cuts for density, temperature and 
magnetic field for W ext = 10% and /3 e ,iono = 0.25/? e . The density profile has a maximum 
value of about n « 7.2 with a characteristic width of the peak of 3 x ri 0 , that is higher 
than the observational value, n as 5.1. The temperature profile shows a two-peak 
distribution. The maximum values in the peaks are T/T 0 rs 3 and T /Tq rs 2.1 that 
correspond to the observational data (Frank et al., 1996; Kivelson et ah, 1996). The 
distance between peaks is about 3 x ri 0 that is a little bit wider than observed value, 

2 x ri 0 . The magnetic field profile has a minimum value in the plasma wake of B « 0.5 
that also corresponds the minimum value of the B y in the observational data. Figure 16, 
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(dotted line) shows the one-dimensional cuts for density, temperature and magnetic 
field for W ext = 30% and fi e ,i ono = 0.25/3 e - The density profile has a maximum value 
of about n « 5.83 with a characteristic width of the peak of 4 x ri 0 , which is a little 
bit wider than that in Fig. 16, (solid line) and is a little higher than observed. The 
temperature profile also shows a two-peak distribution. The maximum values in the 
peaks are T/To ~ 6.8 and T /T 0 ~ 2.4. Note that the first maximum in the temperature 
is higher than observed (Frank et al., 1996; Kivelson et ah, 1996). The magnetic field 
profile shows a minimum value of the magnetic field in the plasma wake B « 0.56 that 
corresponds the observational data. 

3.4. Effect of Ionospheric Conductivity on the Plasma Environment 

In the previous cases we have assumed a high conductivity for the ionosphere and 
Io’s body. The realistic models of Io may include a conducting core, surrounded by a 
poorly conducting mantle (a Moon-like model), or Io may be considered as a poorly 
conducting body. In this section we shall model the Io as a poorly conducting body. 
Note that all other parameters in this case are the same as that shown in Fig. 16, (solid 
line) except the higher diffusion length. Figure 16 (dashed line) shows one-dimensional 
cuts for the density, temperature and magnetic field for the case with a diffusion length, 
ld,i<mo — 0.025. Note that the lines in Fig. 16 corresponds the simulation for 
different sizes of the computational domain. The value of the peak density is 
about n « 7.5 with the same thickness of the peak, « 3r Io as in Fig. 16, (solid line). 
The maximum values of the peak temperature, 3 and 2 are approximately the same as 
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in Fig. 16, (solid line) (see Fig. 16, (dashed line)). However, the magnetic field profile 
shows much stronger field variation in the external region of the plasma wake than in 
Fig. 16, (solid line). The analysis of the magnetic field inside Io shows that the reduced 
conductivity (Zd,iono = 0.025, Fig.18 (dashed line)) near ionosphere and inside Io may 
reduce the magnetic field inside Io by 20 — 40% and the peak plasma density at the 
Galileo trajectory by 5% in comparison with the model with Zd.iono = 0.0025). 

3.5. Comparison with Galileo Observational Data 

The result of the measurements by the particle and field instruments on the Galileo 
Orbiter during the December 1995 flyby of Io provided new and important information 
with which realistic simulations for the plasma interaction can be tested. Along that 
trajectory physical signatures of the wake were seen as a broad depression in the 
magnetic field (Kivelson et al., 1996), sharp peaks in the ion (Frank et al., 1996) and 
electron (Gurnett et al., 1996) densities, a slowing of the plasma in the core of the wake, 
a deep ion temperature decrease in the center of the wake, and a large (factor of 3) 
ion temperature rise in the flanks of the wake (Frank et al., 1996). The magnetic field 
perturbation was broader spatially than the density peak, and showed a double-reversed 
structure, whereby the perturbation (as defined by the difference from the outer Jovian 
B-field value) was actually weaker right near the close-approach point than it was 
somewhat adjacent to the center of the wake. For comparison of the our computational 
model with observation data we made a run with the following plasma parameters: 
/3e,pi == 0.125/3 e>U p, /3e,iono — /5e,up> u.j 0no = 100rio, -Hatmos = 0.06ri o , W ex t — 5%, and 
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Z d>up = 0.00125. The conductivity of Io is chosen so that Zd,io = 0.0125. 

Let us consider first a case in the absence of charge exchange processes (case nee 
Table 1). The total ion production rate was Qion = 4.8 x 10 27 s _1 . Fig. 17 shows the 
comparison of plasma parameters and magnetic field components from the simulation 
model and observation along the Galileo trajectory (Pass 10). Note that we used the 
interpolation of the grid values of the plasma parameters and the magnetic 
field into the position of Galileo spacecraft. We can see that the density peak 
is slightly shifted to the right in comparison with observed one. The left peak in the 
temperature profile is a little bit narrower than the observed one while the right peak 
is a little bit higher than the observed one. The simulation yields a little bit smaller 
value of By at large x, x >> r Io , due to perturbation of the electromagnetic field in the 
region above ( y > 0) and below ( y < 0) the equatorial plane. Unlike the MHD models 
the hybrid model yields a B y profile with reverse structure in the middle of the plasma 
wake. This type of behavior may be explained by the diamagnetic and accelerational 
drift currents in the plasma wake, which are modeled naturally in our kinetic description 
for ions. The maximum value of these currents is located near the equatorial plane 
in the boundary layer (interface) between the external plasma flow and pickup ions 
in the plasma wake. The B z profiles in these simulations correspond very well to the 
observed data however the left maximum in the magnetic field is smaller than observed 
one. Unfortunately fluctuations in magnetic field components are not very small due to 
“shot” noise. So the lack in agreement between observed and computed values of B x is 
probably just the result from less than perfect plasma simulation model (“shot” noise), 
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and the simplified models of the ionosphere and Io’s body. 

Let us consider now the cases with a charge exchange rate that corresponds to 
the maximum value of the density of the atmosphere, n atmos = 10 8 cm -3 . We assume 
here that the charge exchange process is due to only the exponential part of the neutral 
atmosphere. 

Fig. 18 (solid line) shows the comparison of plasma parameters and magnetic field 
components from our simulation model and the observations in the case with a total ion 
production rate, Q j on = 2.27 x 10 27 s -1 . In this case the density profile has a smaller 
effective value but the total temperature is a little bit higher. The left peak in the 
temperature profile is also a little bit narrower than in the observational data while the 
right peak is a little bit higher than the observed one. The perturbation in B y is a little 
bit smaller than in the case of absence of charge exchange, (cf. Fig. 18 (solid line) and 
Fig. 17). 

In case of higher total ion production rate, Q\ OD = 2.57 x 10 27 s -1 (Fig. 18, (dotted 
line)) the density profile has a higher effective value but the total temperature is a little 
bit lower. We can see that the density profile has a depletion on the left side and the 
width of the peak is smaller than in the case with no charge exchange, (cf. Fig. 18, 
(dotted line) and Fig. 17). The left peak in the temperature profile is also a little bit 
narrower than in the observational data while the right peak is a little bit higher than 
the observed one. The perturbation in B y is a little bit smaller than in case of absence 
of charge exchange (cf. Fig. 18, (dotted line) and Fig. 17). 

In a case with a higher charge exchange rate that corresponds to the maximum 
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value of the density of atmosphere, n atmos = 5 x 10 8 cm -3 and a total ion production 
rate, Q[ 0n = 2.1 x 10 27 s _1 we have a satisfactory agreement in the total density profile, 
a smaller total temperature and no agreement in the magnetic field component B y . 

Let us consider the cases with a high charge exchange rate that corresponds to the 
maximum value of the density of atmosphere, n atmos = 10 9 cm~ 3 and an ionization rate 
of Qion = (1-4 - 2.1) x 10 27 s~L In case of low ionization rate, Q ion = 1.4 x 10 27 s -1 , 
the density peak is higher and thinner than the observed one whereas the temperature 
profile is in a good agreement with the observation. The magnetic field B y profile is 
only in qualitative agreement with observation. In the case of a higher ionization rate, 
Qion — 2.1 x 10 27 s~ 1 , the density has a strong double peak structure but the temperature 
profile is in agreement with the observations. However, the magnetic field profile has no 
agreement with observational data. 

In previous cases we assumed that the charge exchange process is only due 
to the lower altitude, exponential part of the neutral atmosphere. Let us consider 
now the results of simulation that includes also the charge exchange due to the 
extended atmosphere which is distributed as r~ 2 . Figure 19 shows the profiles of 
the total density, temperature and magnetic field for different maximum values of 
the density of atmosphere, n atm os = (5 x 10 7 — 10 9 )cm" 3 , and an ionization rate, 

Qion = (3.03- 4.04) x 10 27 s~L 

Let us consider first the cases with a low maximum value of the density of 
atmosphere, n atm os = 5xl0 7 cm“ 3 , and an ionization rate ofQ ion = (3.03-4.04) xlO 27 s' 1 . 
In case of lower ionization rate, Q ion = 3.03 x 10 27 s~ 3 , the density profile is in agreement 
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with the observation, but the temperature is a little bit higher in the peaks (Fig. 19, solid 
line). The magnetic field profile ( B y ) is in a qualitative agreement with the observation 
(Fig. 19, solid line). In the case of higher ionization rate, Q ion = 4.04 x 10 27 s _1 , the 
density profile is a little bit higher in peak than the observed one and the temperature 
is a little bit lower in the left peak. The magnetic field profile ( B y ) is much smoother 
than in the observation (not shown in Fig. 19). 

In the cases with a higher maximum atmosphere density, n atmos = 10 8 cm -3 , the 
maximum value of the wake density becomes much higher than in the observation. 
Although the temperature profiles are in good agreement with the observations, the 
magnetic field component B y is much smoother in the case of Q ion = 3.03 x 10 27 s _1 
(Fig. 19, (dotted line)), and it is in qualitative agreement with the observations in the 
case with Q; on = 4.04 x 10 27 s _1 (not shown in Fig. 19). The increase in the atmosphere 
maximum to n atmos = 5x10® s -3 results in an increase of the density to more than 
twice the observed value. In this case we have Qj on = 4.04 x 10 27 s _1 with a total charge 
exchange rate of Q exch = 5.8 x 10 28 s _1 , Fig. 19, (dashed line). 

In case of a high maximum atmosphere density, n atmos = 10 9 cm -3 there is 
no similarity between simulation and observation. We simulated the cases with 
ionization rates, Qion = (2.1 — 4.04) x 10 27 s _1 . In the case of low ionization rate, 

Qion = 2.1 x 10 27 s -1 , the density profile has two thin peaks (Fig. 19, (dot-dashed 
line)) while the magnetic field profile ( B y ) has a strong variation (Fig. 19, (dot-dashed 
line)) that does not correspond the observation. In case of moderate ionization rate, 
Qion = 3.03 x 10 27 s _1 , the density profile has one peak with a value, ntotai > 10 uq 
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that is much higher than observed one whereas the magnetic field profile ( B y ) is much 
smoother than observed one. In the case of a high ionization rate, Qion = 4.2 x 10 27 s _1 , 
the density profile has one peak with a value, n to tai ~ (8 — 10)n 0 that is much higher 
than observed whereas the magnetic field profile (B y ) is much smoother than observed. 

One of the main issues of the plasma torus - Io interaction is the question about 
the real obstacle in this interaction. There are two possible candidates for this problem 
- (1) mass loading by the pickup ions originally produced from Io’s atmosphere by the 
ionization processes or charge exchange processes between the incoming plasma torus 
ions and the neutral atmosphere, or (2) charge exchange processes between the pickup 
ions originally produced from Io’s atmosphere and the neutral atmosphere. Table 
1 gives a summary of results of simulations for different models that are 
considered in this section. 

In the case of the models with a high density of Io’s atmosphere (10 9 cm -3 ) the 
charge exchange rate becomes very high, Qexch ~ (6.7 — 9) x 10 28 s _1 (cases ceie6, 
ceie7 and seie8, Table 1). Note that in these cases the charge exchange rate for 
incoming ions is much smaller than for pickup ions. We compared the charge exchange 
rate in the external (r > 4ff atmos ) and the internal (r < 4 x 77 atmos ) region. The charge 
exchange rate for incoming ions in the external region has approximately the same value 
as a charge exchange rate in the internal region. However, the charge exchange rate for 
pickup ions in the external region is approximately 20 times smaller than the charge 
exchange rate in the internal region. 

The magnetic field profile ( B y ) does not match the observational data very well 
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for a small ionization rate, Qj on = 2.1 x 10 27 s -1 (case ceie3, Table 1), whereas the 
total density value in the plasma wake is much more than observed one in the case of 
higher ionization rate, Q i 0n = (3.03 - 4.04) x 10 27 s~ x (cases ceiel and ceie2, Table 
1). Hence, the models with a high charge exchange rate, Q exC h = (6.5 - 9) x 10 28 s _1 , 
are not realistic because they cannot explain the observational data. 

The model with n atmos = 5 x 10 8 cm -3 (case ceie5, Table 1), gives the total charge 
exchange rate Qexch = 5.8 x 10 28 s _1 . The models with n atm0 s = 10 8 cm -3 (cases ceie3 
and ceie4, Table 1), give total charge exchange rates of Q ex ch = (1 09- 1.35) x 10 28 s _1 
And, finally, the models with n atmos = 5 x 10 7 cm~ 3 , give total charge exchange rates of 
Qexch = (6.12 - 6.54) x 10 27 s _1 (cases ceiel and ceie2, Table 1). 

The analysis of the models shows that two cases that have a good fit for 
observational data. In the model without charge exchange (case nee, Table 1) the 
simulation results describe well enough the main observational data - a total density, a 
total temperature and the magnetic field ( B y ) (Fig. 17). The more realistic model with 
n atmos = 5 x 10 7 cm -3 , and ionization rate Qj on = 3.03 x 10 27 s _1 (case ceiel, Table 1), 
also fits satisfactorily the observational data (Fig. 19, solid line). This model gives the 
charge exchange rates Q p , e xch = 3.16 x 10 27 s -1 and Q c ,exch = 3.13 x 10 27 s _1 for plasma 
torus ions and pickup ions respectively and would correspond to an MHD model with 
total fresh ion mass-loading rate and the charge exchange rate, which contributes to the 
momentum and energy friction terms, of 1.2 x 10 28 s _1 (Combi, Gombosi and Kabin, 
2002). So, we find that the interaction is dominated in roughly equal parts between 
primary torus ion charge exchange and secondary pickup ion charge exchange. 
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4. Conclusions 

3D Boltzmann simulations of the interaction of the Jovian plasma torus with Io, 
have demonstrated several new features: 

• The whistler (lead front) and quasi-stationary Alfven waves (wings) strongly affect 
plasma flow around Io. 

• The effect of the finite ion gyroradius that results in formation of the asymmetrical 
boundary layer in the vicinity of Io’s ionosphere is important. The plasma 
parameters have a strongly asymmetrical distribution across Io’s wake. The 
kinetic behavior of ion dynamics reproduces the inverse structure of the magnetic 
field (due to drift current) which cannot be explained by standard MHD or 
electrodynamic simulations which do not account for anisotropic ion pressure. 
The diamagnetic effect of non-isotropic gyrating pickup ions broadens the 
B-field perturbation and produces increased temperatures in the flanks of the 
wake, as observed by the Galileo spacecraft, but not explained by previous 
simulations. Note that two-fluid simulation (Saur et al., 2002) produces 
the double-peak signature with spatial scale much smaller than in 
observation. 

• The cold, dense wake is produced, as in MHD, but was not produced in 
electromagnetic simulations without ad hoc addition of bidirectional electrons. 


• The values of temperatures of the electrons which are created and cooled by 
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collision with neutrals in the exosphere and inside the ionosphere may strongly 
affect the pickup ion dynamics along the magnetic field and consequently the pick 
up distribution across the wake. In the absence of an observed global picture of the 
plasma distribution, the simulation serves to demonstrate the wide range of global 
configurations that are possible for various electron temperature descriptions. A 
full hybrid simulation for the plasma including an electron temperature treatment 
that accounts for electron-neutral collision and a detailed neutral description 
would be required to produce an accurate global picture from the first principles. 

• The effective conductivity of Io’s ionosphere may change strongly the distribution 
of the magnetic field near and inside Io. The reduced conductivity (i l dif = 0.025) 
near the ionosphere and inside the Io may reduce the magnetic field inside Io by 
20 — 40% and the peak of the plasma density at the Galileo trajectory by 15%. 

• The best models that fit the observational data well are: the model without charge 
exchange and ionization rate Qi 0n = 4.8 x 10 27 s -1 , and the model with ionization 
rate Q ion = 3.03 x 10 27 s _1 and the total charge exchange rate Q e xch = 6.3 x 10 27 s -1 . 
The best MHD model had a total fresh ion mass-loading rate of 6 x 10 27 s -1 and 
a charge exchange rate, which contributes to the momentum and energy friction 
terms, of 1.2 x 10 28 s -1 (Combi, Gombosi and Kabin, 2002). The reason for the 
factor of 2 difference results from the above-described pressure anisotropy effect 
in the kinetic simulation, and the resulting tighter plasma distribution near the 


equator plane defined by Io. 
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• In the presented calculations the value of the ion gyroradius was much 
larger than the realistic value. So we have to investigate the effects of 
realistic values of the ion gyroradius in future simulations. We expect 
that the smaller ion gyroradius may result in smaller asymmetry of 
the global picture. However, we do not expect the strong changes in 
the plasma and the magnetic field profiles because the used value of 
the ion gyroradius is much smaller than the characteristic scale of the 
problem - the radius of Io. These profiles are controlled primary by the 
inter-penetration of the torus plasma and the pickup ions. 
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5. Figure captions 

Figure 1 . Galileo trajectory close to Io and the system of coordinates. 

Figure 2. Incoming (top) and pickup (bottom) ion density in the x-z plane. The 
case with Q ion = 3.03 x 10 27 s _1 , n atmos = 5 x 10 7 cm~ 3 , /3 e<P[ = /3 e , /3 e>iono = 0.25 /3 e , 
W ext = 5%, and H atmos — 0.06. See explanation in Fig. 6. 

Figure 3. Incoming (top) and pickup (bottom) ion density in the y-z plane for the 
same parameters as Fig. 2. See explanation in Fig. 7. 

Figure 4. Electric (top) and magnetic (bottom) fields in the x-z plane for the 
same parameters as Fig. 2. Figure shows a strong asymmetry of the electromagnetic 
field due to finite gyroradius effects. 

Figure 5. Electric (top) and magnetic (bottom) fields in the y-z plane for the same 
parameters as Fig. 2. Figure shows the formation of an Alfven wing in the direction of 
the main magnetic field. 

Figure 6. Incoming (top) and pickup (bottom) ion velocity arrows in the x-z plane 
for the same parameters as Fig. 2. The figure demonstrates a flow of pickup ions from 
the “corona” across the magnetic field. The incoming ions flow around the effective 
obstacle that is produced by pickup ions and ionosphere. 

Figure 7. Incoming (top) and pickup (bottom) ion velocity arrows in the y-z 
plane for the same parameters as Fig. 2. The figure demonstrates a strong expansion of 
pickup ion “corona” along the magnetic field line. The incoming ions flow around the 


region of extended “corona” . 
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Figure 8. Two-dimensional section for total density in the x-y plane for the same 
parameters as Fig. 2. 

Figure 9. Total density, temperature and magnetic field along the x-axis. The 
case with Pe,pi — 0, Pe,i ono ~~ IF ext ~ 0%, and H a tmos — 0.06. 

Figure 10. Two-dimensional section for total density in the x-y plane. The case 
with Pe iPI = 0, Pe^ono = 0, W ext = 0%, and H atmos = 0.06. 

Figure 11. Total density, temperature and magnetic field along the x-axis. The 
case with Pe,pi — Pet Pe,iono — 0 , W ex t — 0%, and H a tmos — 0.06. 

Figure 12. Two-dimensional section for total density in the x-y plane. The case 
with Pe,pi = Pet Pe,iono = 0, W ext = 0%, and H atmos = 0.06. 

Figure 13. Total density, temperature and magnetic field along the x-axis. The 
Case with Pe,PI = Pet Pe,iono = Pet W ex t = 0%, H a tmos ~ 0.06. 

Figure 14. Two-dimensional section for total density in the x-y plane. The case 
with Pe,PI — Pet Pe,iono = Pet ^ext = 0%, Hatmos ~ 0.06. 

Figure 15. Two-dimensional sections for pickup ion density in the x-y plane. 
The cases with (top) p giPI = 0, p e ,imo = 0; (middle) p e<PI = P e , p e , iono — 0; (bottom) 
Pe,pi = Pet Pe.iono = 0.25 P e . W ext = 0%, and H atmos = 0.06. Figure demonstrates 
asymmetry of pickup ion density across the wake. 

Figure 16. Total density, temperature and magnetic field along the x-axis for 
Pe,pi = Pet Pe,i<mo = 0.25/5 e , H atmos = 0.06, l d ,u P = 0.0025. (solid line) W ext = 10%, 
ld,iono = 0.0025; (dotted line) W ext = 30%, ld,iono = 0.0025; (dashed line) W ext — 10%, 
ld,iono — 0.025. 
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Figure 17. Comparison of Galileo PLS (10 pass) [ Frank et al, 1996] and MAG 
[ Kivelson et al., 1996] (the open circles) data with the Io hybrid model results in absence 
of charge exchange. Q = 4.8 x 10 27 s -1 , /3 e ,p/ = 0.125 /3 e , f3 e ,iono = 0.125/? e , W ext = 5%, 
Hatmos = 0 06, and n iono = 800(100)n o . Id, up = 0.00125 and ldj 0 = 0.0125. 

Figure 18. Same as Fig. 17, but with charge exchange (Sex P^~^j), 
n atmos — 10 8 cm -3 : (solid line) Q = 2.27 x 10 27 s _1 ; (dashed line) Q — 2.57 x 10 27 s _1 . 

Figure 19. Same as Fig. 17, but with charge exchange (^ + Sexp ) : (solid 
line) n atmos = 5 x 10 7 cm~ 3 , n iono = n 0 , Q = 3.03 x 10 27 s _1 , Q v , exc h = 3.169 x 10 27 s -1 , 
Qc,exch = 3.13 x 10 27 s -1 , Idjo = 0.00125; (dotted line) n atmos = 10 8 cm“ 3 , 
niono = n 0 , Q = 3.03 x 10 27 s _1 , Q v , exc h = 5-27 x 10 27 s _1 , Q c ,exch ~ 5.6 x 10 27 s -1 , 
h,io = 0.00125; (dashed line) n atmos = 5 x 10 8 cm -3 , ni ono — n 0 , Q = 4.04 x 10 27 s _1 , 
Qp,exch = 1-83 x 10 28 s _1 , Q c ,exch = 3.97 x 10 28 s _1 , Idjo = 0.0125; (dot-dashed 
line) n a t mos = 10 cm , n^ono — no, Q ~ 2.1 x 10 s , Q p ,exch — 2.95 x 10 s , 
Qc,exch = 3.75 x 10 28 s _1 , l d jo = 0.0125. 



Table 1 . Dependence of torus ion charge exchange rate and pickup ion charge exchange rate 
versus the maximum value of atmospherical density and ionization rate (Q 0 = 10 27 s -1 ). 


* denote the regimes that fit satisfactory the observational data 

CaSe lQ^cm- 3 Qp 4o Ch Qc $o Ch Comments 

model without charge exchange 

nee 4.8 0 0 a density is a little bit wide and it is in good agreement with an observation, 

temperature and magnetic field are in a good agreement with an observation, Fig. 17* 
models with charge exchange in the internal region (neutral ~ B exp Hfl ) 

“atmos 

ceil 1 2.27 density, temperature and B y in satisfactory agreement with observation, Fig. 18, 

(solid line) 

cei2 1 2.57 density and temperature in satisfactory agreement with observation, B y component 

is a weaker than observed, Fig. 18, (dashed line) 

cei3 5 2.1 density in good agreement with observation; temperature and B y are not in agreement 

with observation 

cei4 10 1.4 12.6 14 density with strong narrow peak, temperature in satisfactory agreement with 

observation; B y is weaker than observed 

cei5 10 2.1 12 20.2 density with two peaks, temperature is in satisfactory agreement with observation; 

By has much stronger depletion than in observation 

models with charge exchange in the internal and external region (n neutra i k, 4 + Bex p r u ~ r k » ) 

r ** atm os 

ceiel 0.5 3.03 3.169 3.13 density, temperature and B y in good agreement with observation. Fig. 19*, (solid line) 

ceie2 0.5 4.04 2.67 3.87 density and temperature in satisfactory agreement with observation, B y is a weaker 

than observed 

ceie3 1 3.03 5.27 5.6 density is higher than the observation, variation in B y is a weaker than observed, 

Fig. 19, (dotted line) 

ceie4 1 4.04 5.55 7.93 density is higher than observation, variation in B y is weaker than observed 

ceie5 5 4.04 18.3 39.7 density is much higher than in an observation, variation in the B y is strong but profile 

is wider than observed, Fig. 19, (dashed line) 

ceie6 10 2.1 29.5 37.5 two-peak density profile, variation in the B y is much stronger than observed, 

Fig. 19, (dot-dashed line) 

ceie7 10 3.03 28. 44.8 density is much higher than in observation, variation in the B y is much weaker than 

observed 

ceie8 10 4.04 27.5 60. density is much higher than in observation, variation in B y is much weaker than 


observed 
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Figure 1 . Galileo trajectory close to Io and the system of coordinates. 





Figure 2. Incoming (top) and pickup (bottom) ion density in the x-z plane. The case 
ceiel (Table 1) with /9 e .pi = A>,iono = 0.25 /3 e , W ext = 5%, and I/ atmos = 0.06. See 
explanation in Fig. 6. 
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Figure 3. Incoming (top) and pickup (bottom) ion density in the y-z plane for the same 
parameters as in Fig. 2. See explanation in Fig. 7. 
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Figure 4. Electric (top) and magnetic (bottom) fields in the x-z plane for the same 
parameters as in Fig. 2. Figure shows a strong asymmetry of the electromagnetic field 
due to finite gyroradius effects. 








Figure 6. Incoming (top) and pickup (bottom) ion velocity arrows in the x-z plane for 
the same parameters as in Fig. 2. The figure demonstrates a flow of pickup ions from the 
“corona” across the magnetic field. The incoming ions flow around the effective obstacle 
that is produced by pickup ions and ionosphere. 





Figure 7. Incoming (top) and pickup (bottom) ion velocity arrows in the y-z plane for 
the same parameters as in Fig. 2. The figure demonstrates a strong expansion of pickup 
ion “corona” along the magnetic field line. The incoming ions flow around the region of 


extended “corona”. 
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Figure 8. Two-dimensional cross section for total density in the x-y plane at z = 1.5rj 0 
for the same parameters as in Fig. 2. 
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Figure 9. Total density, temperature and magnetic field along the x-axis at z = 1.5rj 0 . 
The case with /? ejP i = 0, /3 e ,iono = 0, W ext = 0%, and i^atmos = 0.06. 
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Figure 10. Two-dimensional cross section for total density in the x-y plane at z = 1.5ri 0 . 
The case with /3 e>P i = 0, /3 e , i ono = 0, W ext = 0%, and // atm0 s = 0.06. 
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Figure 11. Total density, temperature and magnetic field along the x-axis at z = 1.5ri 0 . 
The case with &,pi = ft, ft,k>no = 0, W^t = 0%, and tf atmos = 0.06. 
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Figure 12. Two-dimensional cross section for total density in the x-y plane at z = 1.5ri 0 . 
The case with /J e ,Pi = ft, ft,iono = 0, W ext = 0%, and // a tmos = 0.06. 
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Figure 13. Total density, temperature and magnetic field along the rr-axis at z — 1.5ri 0 . 
The case with y3 e ,pi = Pe, Pe,iono = Pe, W ext = 0%, # a tmos = 0.06. 
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Figure 14. Two-dimensional cross section for total density in the x-y plane at z = 1.5ri 0 . 
The case with /3 e ,pi = Pe, Pe,\ono = Pe, W ex t = 0%, H itmos = 0.06. 






Figure 15. Two-dimensional cross sections for pickup ion density in the x-y plane at 
2 = 1.5r Io . The cases with (top) &,pi = 0, /? e ,iono = 0; (middle) &.pi = /3 e ,iono = 0; 
(bottom) /? e ,pi = p e , Pejono = 0.25/? e . W ext = 0%, and ifatmos = 0.06. Figure demon- 


strates asymmetry of pickup ion density across the wake. 
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Figure 16. Total density, temperature and magnetic field along the x- axis at, z = 1.5r Io 
for cases with /? e ,Pi = Pe, Pe,iono = 0.25&, H & tm os = 0.06, / djUp = 0.0025. (solid line) 
W ext = 10%, Zd.iono = 0.0025; (dotted line) W ext = 30%, /d.iono = 0.0025; (dashed line) 
W ext = 10%, l d} iono = 0.025. 
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Figure 17. Comparison of Galileo PLS (10 pass) [Frank et al., 1996] and MAG [Kivel- 
son et al., 1996] (the open circles) data with the Io hybrid model results along the 
trajectory in the absence of charge exchange (case nee, Table 1). /3 e ,pi = 0.125/3 e , 
fiejono = 0.125 fi e , W e x t = 5%, # a tmos = 0.06, and niono = 800(100)n 0 . id, up = 0.00125 and 




^d,io = 0.0125. 
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. 17, but with charge exchange ( B exp hJ 1 I s ), n atmos = 10 8 cm 2 3 * : 

>le 1); (dashed line, case cei2, Table 1). 
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Figure 19. Same as Fig. 17, but with charge exchange {£+B exp ): (solid line, case 
ceiel, Table 1) nj 0n0 = no, id,io = 0.00125; (dotted line, case ceie3, Table 1) /d,i 0 = 0.00125; 
(dashed line, case ceie5, Table 1) Zd,io = 0.0125; (dot-dashed line, case ceie6, Table 1) 


^d,io = 0.0125. 



